% function fig1(filenumber)
% count=filenumber;
% 
%     N=readdata(count,1);
%     S=readdata(count,2);
%     T=readdata(count,3);
% %     multigridlevels=readdata(count,4);
% %     bound=readdata(count,5);
%     time=readdata(count,6);
% %     dtnom=readdata(count,7);
%     dx=readdata(count,8);
%     dz=readdata(count,9);
%     mx=readdata(count,10);
%     mz=readdata(count,11);
%     
%     
% x=zeros(1,S); z=zeros(1,T);
% for s=1:S-1 x(s+1)=x(s)+dx(s); end
% for t=1:T-1 z(t+1)=z(t)+dz(t); end
    

    
% mTemp=readdata(count,15);
% [grid,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,2,2); ind1=find(grid<1);
% avtemp=0*zgrid;
% for t=1:length(avtemp);
%     avtemp(t)=mean(grid(:,t));
% end
% load colormapfig1a;
% figure(1); clf;
% axes('pos',[0.04 0.54 0.42 0.42]);
% [Zgrid,Xgrid] = meshgrid(zgrid,xgrid); grid(ind1)=nan; pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colormap(cmapfig1a); colorbar; axis([0 2000 -40 670]); title(['Temperature (^{\circ}C)']); xlabel('Distance (km)'); ylabel('Depth (km)'); clear('mTemp');
% 
% axes('pos',[0.04 0.05 0.42 0.42]);
% plot(avtemp,(zgrid-40e3)/1e3); hold on; set(gca,'ydir','rev');
% plot(adiabat(zgrid,40e3),(zgrid-40e3)/1e3,'r'); axis([0 1900 0 670]); xlabel('Temperature(^{\circ}C)'); ylabel('Depth (km)');
% legend('Horisontally averaged temperature ','Theoretical isentropic temperature at potential \newline temperature T_p=1315 ^{\circ}C','location','southwest');
% 

load surface_time;
S=length(x);
timehistint=linspace(timehist(1),timehist(end),length(timehist)); dt=timehistint(2)-timehistint(1);

lowcut=0;
surfacehisthf=surfacehist;
d=5;
qhist=2.5*(temphist(:,:,d+1)-temphist(:,:,d))/(depthrange(d+1)-depthrange(d)); meanq=mean(qhist(:))+(0.5*(depthrange(d+1)+depthrange(d))-40e3)*0.75e-6;
for s=1:S
    tempsurface=interp1(timehist,surfacehist(s,:),timehistint); tempsurface(1)=tempsurface(2);
    temptemp=interp1(timehist,avtemphist(s,:),timehistint);
    tempq=interp1(timehist,qhist(s,:),timehistint);
    avtemphist(s,:)=bandpass(temptemp,lowcut,inf,dt);
    qhist(s,:)=bandpass(tempq,0,inf,dt);
    surfacehist(s,:)=bandpass(tempsurface,lowcut,inf,dt);
    surfacehisthf(s,:)=bandpass(tempsurface,1/50,inf,dt);
    
    for d=1:length(depthrange)
        temptemp=interp1(timehist,temphist(s,:,d),timehistint);
        temptemp=bandpass(temptemp,lowcut,inf,dt);
        temphist(s,:,d)=temptemp;
    end
end
qhist=qhist+meanq;

[TIMEHIST X]=meshgrid(timehistint,x);

subplot(4,4,[1 2 5 6]);
pcolor(X/1e3,TIMEHIST,-surfacehist); shading flat; xlabel('Distance (km)'); ylabel('Time (Ma)'); colorbar; title('Variations in vertical surface position (m)');

subplot(4,4,[3 4 7 8]);
pcolor(X/1e3,TIMEHIST,qhist); shading flat; xlabel('Distance (km)'); ylabel('Time (Ma)'); colorbar; title('Surface heat flux (W/m^2)');

subplot(4,4,[9 10 13 14]);
pcolor(X/1e3,TIMEHIST,-surfacehisthf); shading flat; xlabel('Distance (km)'); ylabel('Time (Ma)'); colorbar; title('Variations in vertical surface position with periods less than 50 Ma (m)');

surfacehistflip=[surfacehist; flipud(surfacehist)];
correlation=zeros(1,S); weight=zeros(1,S);
for refindex=1:S
    for s=refindex:refindex+S-1
        correlation(s-refindex+1)=correlation(s-refindex+1)+(surfacehistflip(s,:)-mean(surfacehistflip(s,:)))*(surfacehistflip(refindex,:)-mean(surfacehistflip(refindex,:)))';
        weight(s-refindex+1)=weight(s-refindex+1)+1;
    end
end
correlation=correlation./weight;
correlation=correlation/correlation(1);
subplot(4,4,[11 12]); plot(x/1e3,correlation); axis([0 2000 -0.5 1]); grid on; title('Normalized lateral correlation of surface variations'); xlabel('Distance (km)');

surfacehistflip=[surfacehisthf; flipud(surfacehisthf)];
correlation=zeros(1,S); weight=zeros(1,S);
for refindex=1:S
    for s=refindex:refindex+S-1
        correlation(s-refindex+1)=correlation(s-refindex+1)+(surfacehistflip(s,:)-mean(surfacehistflip(s,:)))*(surfacehistflip(refindex,:)-mean(surfacehistflip(refindex,:)))';
        weight(s-refindex+1)=weight(s-refindex+1)+1;
    end
end
correlation=correlation./weight;
correlation=correlation/correlation(1);
subplot(4,4,[15 16]); plot(x/1e3,correlation); axis([0 2000 -0.5 1]); grid on; title('Normalized lateral correlation of surface variations with periods less than 50 Ma'); xlabel('Distance (km)'); 


figure; pcolor(X/1e3,TIMEHIST,-(surfacehist-surfacehisthf)); shading flat; xlabel('Distance (km)'); ylabel('Time (Ma)'); colorbar; title('Variations in vertical surface position with periods higher than 50 Ma (m)');